Pulsar recoil by large-scale anisotropies in supernova explosions 
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Assuming that the neutrino luminosity from the neutron star core is sufficiently high to drive su- 
pernova explosions by the neutrino- heating mechanism, we show that low-mode (/ = 1, 2) convection 
can develop from random seed perturbations behind the shock. A slow onset of the explosion is cru- 
cial, requiring the core luminosity to vary slowly with time, in contrast to the burst-like exponential 
decay assumed in previous work. Gravitational and hydrodynamic forces by the globally asymmetric 
supernova ejecta were found to accelerate the remnant neutron star on a timescale of more than a 
second to velocities above 500 kms~^, in agreement with observed pulsar proper motions. 

PACS numbers: 97.60.Bw, 97.60.Gb, 95.30.Jx, 95.30.Lz 



Young pulsars are observed to have average space ve- 
locities of 200-500 kms~^ with highest values above 
1000 kms~^ and still ambiguous hints for a double peak 
distribution H. There is no clear statistical correlation 
vifith the magnetic moment or rotation of the pulsar, al- 
though in two cases (Vela and Crab) the direction of 
motion appears to be aligned with the spin axis. 

A connection of the pulsar motions with the super- 
nova (SN) explosion is suggested by neutron star (NS) - 
SN remnant associations and by the properties of binary 
systems with one or both components being a NS. Natal 
kicks are required, e.g., by the spin-orbit misalignment 
and high orbital eccentricities observed in some binaries 
(for a review, see 0)- 

Various mechanisms have been proposed to explain 
these kicks. One possibility invokes asymmetric mass 
ejection during the SN j^. This may be caused by large- 
scale density inhomogeneities in the pre-collapse core of 
the progenitor star '411 or convective instabilities in the 
neutrino- heated layer behind the SN shock ^ ^ . The 
kicks could also be a consequence of unequal momentum 
fluxes in a jet and an anti-jet that might be linked to the 
start of the SN explosion or to the NS formation 

Alternatively or in addition, anisotropic neutrino (v) 
emission from the nascent ("proto-") neutron star (PNS) 
might transfer the momentum [^. A "neutrino rocket 
engine" of the latter kind could result from the magnetic 
field-strength dependence of i/-matter interactions if ex- 
tremely strongfields with hemispheric asymmetries build 
up in a PNS |9(. The v transport also depends on the 
field direction, e.g., through parity violating corrections 
of weak interaction cross sections [liJ] or due to resonant 
flavor transitions ■ In case of a significant dipole com- 
ponent such effects can lead to kicks of a few 100 kms~^ 
for magnetic fields in excess of 10^^ G Q. 

In this Letter we present new two-dimensional (2D) 
calculations of hydrodynamic instabilities during the on- 
set of SN explosions which show that global asymmetries 



and the PNS recoil can naturally grow to a sufficient 
size without invoking artificial initial conditions, extreme 
physical assumptions, or exotic v physics. Our computa- 
tions improve previous ones |^ with respect to numerical 
resolution, a full 180° lateral grid, and the treatment of 
v transport, extending them also in the computed evolu- 
tion time and model set. 

Modeling concepts. We assume that the explosion is 
powered by rz-energy deposition between the PNS and 
the SN shock T^l . Although the currently most elaborate 
numerical models arc not able to confirm the viability of 
this I'-hcating mechanism, it is still the best-studied and 
most promising way to explode massive stars |13]. 

SN theory is currently hampered by our incomplete 
knowledge of the nuclear physics and u interactions in the 
dense matter inside the PNS. This implies uncertainties 
for self-consistent models of the full problem, e.g. with 
respect to the magnitude of the i/ fluxes emitted by the 
cooling PNS. Therefore, we replace the shrinking, high- 
density core of the PNS by a gravitating sphere whose 
radius coincides with the contracting, impenetrable inner 
boundary of our computational grid. 

At this boundary, number and energy fluxes of v and 
D of all three lepton flavors are imposed with chosen ini- 
tial values and time dependence. In all simulations these 
boundary values were taken to be isotropic and were kept 
constant for a chosen period of time after shock forma- 
tion (either 0.7s or Is, within which 1/3 of the gravi- 
tational binding energy of the nascent NS was assumed 
to be radiated away in neutrinos). The grid boundary is 
located somewhat below the neutrinosphere of Ve- It has 
an optical depth which, for typical v energies, rises from 
few initially to several 100 at the end of the simulations. 
The use of this inner boundary allows us to explore the 
response of the collapsing SN core to different v lumi- 
nosities. Higher values of the latter lead to larger energy 
deposition behind the SN shock and therefore to more 
powerful explosions. The dominant heating reactions are 
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FIG. 1: Explosion timescale tcxp and, at one second after SN 
shock formation, explosion energy i5exp, NS baryonic mass 
Mns, gas anisotropy parameter agas, NS recoil velocity Vns, 
and NS acceleration Uns (from bottom to top) vs Ue plus i>e 
luminosity. Lib, at the inner boundary. The symbols corre- 
spond to different models of 15 Mq stars, the triangles to a 
case including rotation (see text for details). 



Ue and Pe absorption on free nucleons. 

The "lightbulb" approximation employed previ- 
ously ignored time retardation effects and did not 
take into account radial variations of the fluxes due to 
z/-matter interactions in the cooling and heating lay- 
ers between PNS and SN shock. To improve on that 
in the present work, we make use of the zeroth mo- 
ment of the Boltzmann transport equation in the form 



dtL + cdrL = 47rr^cQ~ — kcL . Here L = L{r,t) is 
the total v number flux or j/ luminosity, respectively, 
the rate of v loss by the stellar medium per unit vol- 
ume, K = Anr^ c / (Lc) the corresponding absorptiv- 
ity, c the speed of hght, and c the "effective speed" of v 
propagation, which is governed by diffusion at high den- 
sities and reaches c at large radii. The integration for L 
in radius r and time t (lateral flux components are ig- 
nored) can be done analytically when Q~ , k and c are 
assumed to be constant within the cells of the numer- 
ical grid (consistent with this, the above equation was 
derived by employing dtc = 0). Instead of determining 
c from the solution of the Boltzmann equation, we use 
an analytic representation in terms of the optical depth 
that was obtained from fitting results of detailed v trans- 
port in the outer layers of the SN core. Neutrino-matter 
interactions via charged-current processes with nucleons, 
thermal pair creation, and scattering off nuclei, n, p, e~, 
and e"*" are evaluated by adopting Fermi-Dirac i' spectra 
with a temperature determined by the ratio between i' 
energy and number density. The chemical potentials of 
the spectra are taken to be equal to the equilibrium val- 
ues at high optical depths and approach values near zero 
outside of the neutrinospheres. 

This approximate v transport retains the hyperbolic 
character of the transport problem, conserves lepton 
number and energy globally, and reproduces basic prop- 
erties of accurate Boltzmann transport calculations de- 
spite of radical simplifications. It is coupled to our 2D hy- 
drodynamics code by operator-splitting with a predictor- 
corrector step. 

Our hydrodynamics code and equation of state are de- 
scribed in Ref. {l4]- We typically use 400 geometrically 
spaced radial zones and one degree lateral resolution. 
The V transport is solved in each angular bin separately. 
The 2D simulations were started some milliseconds after 
SN shock formation from detailed core-collapse models 
with fluid velocities randomly perturbed with an ampli- 
tude of typically 0.1%. 

Results. Figure ^ shows four sequences of runs that 
followed the post-bounce evolution of different 15 Mq 
progenitors for one second. The crosses mark results 
based on the use of Model WPE15 0, the solid dots of 
Model LSC15 [13, the open circles of Model sl5s7b2 [l^, 
and the triangles of a model with a structure like the lat- 
ter but including rotation Starting with a rotation 
period of 12 s in the iron core 0|, the post-collapse core 
spins differentially within 10-20 ms and speeds up as it 
contracts. The rotation period in the j/-heating layer 
varies between some 10ms and several 100ms |13l |. 

Model LSC15, in particular, differs from the others by 
having significantly higher densities at the edge of the 
iron core and in the silicon shell (at a time when the 
cores have evolved to the same central density). This 
delays the start time, texp, of the explosion, reduces the 
explosion energy, E^xp , and leads to a larger NS baryonic 
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FIG. 2: Two mod- 
els (based on progeni- 
tor WPE15) at Is af- 
ter SN shock formation. 
The PNS is at the co- 
ordinate center. The 
left plot (scale reduced 
to show more details) 
displays a case with 
Lib = 2.97x10^^^ ergs-\ 
Scxp = 0.4 X lO'-i erg, 
and a recoil velocity of 
?;ns = -350 kms"^ The 
model on the right has 
Lib =4.45x10^^ ergs-\ 
Eexp = 1.2 X 10^1 erg, 
and VnB = +520 kms~^. 



mass, Mns, for a given value of the boundary luminosity, 
Lib, of i/e plus Pe (Fig. Q). L^cxp IS defined as the total 
energy (internal plus kinetic plus gravitational) of the 
SN ejecta, integrated over all matter where the sum of 
the corresponding specific energies is positive, iexp is the 
post-bounce time when Ecxp reaches 10^^ erg, and AIns 
is the gas mass with densities above 10^^ gcm~'^ plus the 
central point mass at one second. 
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FIG. 3: Explosion energies (bottom) and NS velocities vs time 
for some simulations, showing large acceleration for cases with 
high «ns (thick lines) even at 1 s after bounce. The symbols 
and line styles refer to the different model sequences of Fig. 

Figure ^ reveals the generic trend that a higher lu- 
minosity Lib from the NS core causes the explosion to 
develop faster and to become more energetic. Because 
the period of mass accretion by the PNS is reduced, 
this implies a smaller NS mass. The time until the re- 
vived bounce-shock reaches a certain radius (correlated 
with tcxp) depends sensitively on the progenitor struc- 
ture and the core v luminosity. Rotation systematically 
increases the explosion energy by 20-50%, because cen- 
trifugal forces delay matter from being accreted onto the 



PNS, and thus keep it in the r^- heating region to accu- 
mulate energy by an d absorption. This is basically 
in agreement with Ref. [l3| , where rotation was found to 
stabilize the standing accretion shock at a larger radius. 

The layer between PNS and SN shock is convectively 
unstable according to the Ledoux-criterion because of a 
negative entropy gradient established by v-eneigy depo- 
sition. Within ~50 ms after shock formation, convection 
sets in, supporting the start of the explosion @ . Ini- 
tially the convective cells are small, but they begin to 
merge to larger entities. Three-dimensional (3D) simu- 
lations agree with this 2D result In case of rapid 
shock acceleration convection freezes, and small struc- 
tures characterize the flow pattern until late times. How- 
ever, when the stagnant shock expands slowly, small cells 
have time to merge to very large buoyant bubbles, sepa- 
rated by only a few narrow, supersonic downflows which 
carry low-entropy matter from the shock to the PNS sur- 
face. Moreover, global pulsations can develop with a 
dominance of low-order {I = 1,2) modes. Consequently, 
the density distribution becomes highly aspherical and 
the explosion breaks out with a very large hemispheric or 
polar-to-equatorial asymmetry and corresponding shock 
deformation. In some cases a single long-lasting accre- 
tion funnel was found to persist for a second or even 
longer (Fig. ^ . We emphasize that such structures did 
not preferentially occur along the polar (z-) direction of 
our spherical grid (where a coordinate singularity exists), 
but developed in arbitrary orientations. 

The development of a stable, volume-filling I — 1 mode 
was proposed before |2fll | . It is supported by analytic ar- 
guments for thermal instabilities in fluid spheres [2l| and 
is also observed in 3D hydrodynamic simulations of pul- 
sating, convective red giant stars Determining the 
duration of such a phenomenon in the time-dependent en- 
vironment of an exploding SN requires numerical mod- 
eling. The coherent, low-order oscillations of the fluid 
beneath the shock in our simulations look similar to the 
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recently discovered non-radial instabilities in adiabatic 
flow behind standing accretion shocks j23i| , which can be 
understood in terms of a "vortical-acoustic cycle" 
Doing numerical experiments we found that this phe- 
nomenon might indeed play a role, although Ledoux- 
instability due to i^-heating clearly starts the convective 
activity, and z/-cooling around the neutrinosphere seems 
to damp the energetic amplification of the feedback cycle 
between turbulence and pressure waves. 

Anisotropic mass ejection can be associated with a 
high linear momentum taken up by the compact rem- 
nant. The corresponding asymmetry of the explosion 



is expressed by the parameter 



|-Pz,gas|/-Pg: 



I / dmvz\/ J dm |i7| where the integrals are performed 
over the ejecta mass and Pz.gas is the gas momentum 
along the z-direction of the 2D grid. We find values 
for Qfgas up to 0.33 (Fig. QJ. The NS recoil velocities, 
^'ns = Ogas-Pgas/A^ns, Can be close to zero but also more 
than 500 kms^^. In some cases the acceleration, ttnsj 
continues on a high level beyond the 1 s of computed 
evolution (Figs. ^ ^ and significantly larger terminal 
velocities can be expected. It is mediated by the long- 
range gravitational force between the asymmetrically dis- 
tributed ejecta and the PNS. Hydrodynamic forces play 
a role only as long as downflows reach the PNS, and 
anisotropic v emission contributes insignificantly. There 
is no correlation of Vns with the progenitor model. Rota- 
tion does not inhibit large kicks. The final recoil velocity 
depends stochastically on the initial seed perturbation 
and the highly nonlinear growth of the convective struc- 
tures. There is also no obvious correlation with the ex- 
plosion energy. Fig. |21 (by the thin and thick solid lines) 
and Fig. ^ (coinciding points) demonstrate that nearly 
the same i?exp can be associated with large or small Vns- 
Conclusions. We have shown that globally anisotropic 
mass ejection and NS acceleration can result from convec- 
tive overturn and low-order oscillations of the i^-heated 
layer in a SN core. Low-mode convection turned out to 
develop from random seed perturbations in case of a slow 
onset of the explosion which gives the convective struc- 
tures time to merge. This was disfavored by our previ- 
ous choice of a strongly time-dependent, exponentially 
decaying core v luminosity in Refs. 0Q| but is possible 
with the less burst-like (because constant) Lib assumed 
in this work. Our models suggest a consistent picture in 
which z^-energy deposition can be responsible for the SN 
explosion, for pulsar kicks, and for global asymmetries 
observed in many supernovae. The simulations need to 
be continued to later times to allow for quantitative con- 
clusions on the morphological properties of the ejecta as, 
e.g., inferred from polarization measurements. Statisti- 
cal information about the distribution of intrinsic pulsar 
velocities requires 3D simulations, a better fundamental 
understanding of the explosion mechanism, and a large 
sample of simulations for progenitor stars with different 
masses and rotation rates. The discussed kick mecha- 



nism does not enforce a strict alignment of pulsar spin 
and space velocity, although the rotation axis of a star 
defines a naturally preferred direction, which might dis- 
favor large misalignments. 

We thank S.W. Bruenn and M. Rampp for post-bounce 
core-collapse models and M. Limongi and S. Woosley 
for their progenitor models. Support by the Sonder- 
forschungsbereich 375 on "Astroparticle Physics" of the 
Deutsche Forschungsgemeinschaft is acknowledged. The 
simulations were done at the Rechenzentrum Garching 
and the Interdisciplinary Centre for Computational Mod- 
elling in Warsaw. 



[1] A.G. Lyne and D.R. Lorimer, Nature (London) 369, 
127 (1994); J.M. Cordes and D.F. Chernoff, Astro- 
phys. J. 505, 315 (1998); Z. Arzoumanian, D.F. Cher- 
noff, and J.M. Cordes, Astrophys. J. 568, 289 (2002); 

B. M.S. Hansen and E.S. Phinney, Mon. Not. R. Astron. 
Soc. 291, 569 (1997); C. Fryer, A. Burrows, and W. Benz, 
Astrophys. J. 496, 333 (1998) 

[2] D. Lai, D.F. Chernoff, and J.M. Cordes, Astrophys. J. 
549, 1111 (2001) 

I.S. Shklovskii, Sov. Astron. 13, 562 (1970) 
A. Burrows and J. Hayes, Phys. Rev. Lett. 76, 352 
(1996); D. Lai and P. Goldreich, Astrophys. J. 535, 402 
(2000) 

M. Herant, W. Benz, W.R. Hix, C.L. Fryer, and S.A. 
Colgate, Astrophys. J. 435, 339 (1994); A. Burrows, J. 
Hayes, and B.A. Fryxell, Astrophys. J. 450, 830 (1995) 
H.-T. Janka and E. Miiller, Astron. Astrophys. 290, 496 
(1994); 306, 167 (1996) 

R. Cen, Astrophys. J. 507, L131 (1998); A.M. Khokhlov 
et al, Astrophys. J. Lett. 524, L107 (1999) 
N.N. Chugai, Sov. Astron. Lett. 10, 87 (1984); A. Bur- 
rows and S.E. Woosley, Astrophys. J. 308, 680 (1986) 

G. S. Bisnovatyi-Kogan, Astron. Astrophys. Trans. 3, 287 
(1993); D. Lai and Y.-Z. Qian, Astrophys. J. 505, 844 
(1998) 

C. J. Horowitz and G. Li, Phys. Rev. Lett. 80, 3694 

(1998) ; P. Arras and D. Lai, Astrophys. J. 519, 745 

(1999) ; Phys. Rev. D 60, 043001 (1999) 
A. Kusenko and G. Segre, Phys. Rev. Lett. 77, 4872 
(1996); E. Nardi and J.I. Zuluaga, Astrophys. J. 549, 
1076 (2001) 

H. A. Bethe and J.R. Wilson, Astrophys. J. 295, 14 
(1985) 

R. Buras, M. Rampp, H.-T. Janka, and K. Kifonidis, 
Phys. Rev. Lett. 90, 241101 (2003) 

K. Kifonidis, T. Plewa, H.-T. Janka, and E. Miiller, As- 
tron. Astrophys. 408, 621 (2003) 

S.E. Woosley, P.A. Pinto, and L. Ensman, Astrophys. 
J. 324, 466 (1988); S.W. Bruenn, in Nuclear Physics m 
the Universe, edited by M.W. Guidry and M.R. Strayer 
(TOP, Bristol, 1993), p. 31 

M. Limongi, O. Straniero, and A. Chieffi, Astrophys. J. 
Suppl. 129, 625 (2000) 

S.E. Woosley and T.A. Weaver, Astrophys. J. Suppl. 101, 
181 (1995) 

A. Heger, S.E. Woosley, N. Langer, and H.C. Spruit, 



[6: 
[7: 
[s: 

[9 
[10 

[11 

[12 
[13 
[14 
[15 



[16 
[17 
[18 



5 



astro-ph/0301374 
[19] C.L. Fryer and M.S. Warren, Astrophys. J. Lett. 574, 
L65 (2002) 

[20] M. Herant, Phys. Rep. 256, 117 (1995); C. Thompson, 
Astrophys. J. 534, 915 (2000) 

[21] S. Chandrasekhar, Hydrodynamic and Hydromagnetic In- 
stabilities (Dover, New York, 1981), p. 234 



[22] P.R. Woodward, D.H. Porter, and M. Jacobs, in 3D Stel- 
lar Evolution, ASP Conf. Series, Vol. 293 (ASP, San 
Francisco, 2003), p. 45 

[23] J.M. Blondin, A. Mezzacappa, and C. DeMarino, Astro- 
phys. J. 584, 971 (2003) 

[24] T. Foglizzo, Astron. Astrophys. 392, 353 (2002) 



